####### MARRIAGE ANALYSIS ###############
st_fem = dn.dta %>% filter(st==1 & fem==1) # ST women 

######## TABLE J.20 #####
controls_dn = c("age + sc + st + fem")
base_dn = paste("marry_other_caste_bin ~ vwreserv*vstreserv")
base_dn_scst = paste("marry_other_caste_bin ~ vwreserv*vscstreserv")

dn.lm1 = paste(base_dn) %>% as.formula() %>%
  lm_robust(data = dn.dta %>%
              rename(vwreserv = reserved_fem, vstreserv = reserved_st),
            clusters = uniq_panch_code, se_type = "stata")

dn.lm2 = paste(base_dn) %>% as.formula() %>%
  lm_robust(data = dn.dta %>%
              rename(vwreserv = reserved_fem, vstreserv = reserved_st),
            clusters = uniq_panch_code, se_type = "stata",
            fixed_effects = ~uniq_dist_code)

dn.lm3 = paste(base_dn, "+", controls_dn, "+", "prop_st") %>% as.formula() %>%
  lm_robust(data = dn.dta %>% rename(vwreserv = reserved_fem, vstreserv = reserved_st),
            clusters = uniq_panch_code, se_type = "stata",
            fixed_effects = ~uniq_dist_code)


dn.lm5 = paste(base_dn_scst) %>% as.formula() %>%
  lm_robust(data = dn.dta %>% rename(vwreserv = reserved_fem, vscstreserv = reserved_scst),
            clusters = uniq_panch_code, se_type = "stata")

dn.lm6 = paste(base_dn_scst, "+", controls_dn) %>% as.formula() %>%
  lm_robust(data = dn.dta %>% rename(vwreserv = reserved_fem, vscstreserv = reserved_scst),
            clusters = uniq_panch_code, se_type = "stata",fixed_effects = ~uniq_dist_code)

dv.means = rep(round(mean(dn.dta$marry_other_caste_bin, na.rm = T), 3), 5)

screenreg(l = list(dn.lm1, dn.lm2, dn.lm3, dn.lm5, dn.lm6),
       include.ci = F, include.rmse = F,
       include.nobs = T,include.adjrs = FALSE, include.rsquared = F,
       digits = 3, include.nclust = F,
       omit.coef = c("age|^sc$|^st$|fem|Intercept|prop\\_st"),
       custom.coef.names = c("Women's Quota", "Caste Quota (ST)", "W X C (ST) Quota",
                             "Caste Quota (ST/SC)", "W X Caste (ST/SC) Quota"),
       stars = c(0.01, 0.05, 0.1),
       custom.model.names = c("All", "All", "All", "All", "All"),
       custom.gof.rows = list(
         "Caste Quota Type" = c("ST", "ST", "ST", "ST/SC", "ST/SC"),
         "DV Mean" = dv.means, 
         "Controls" = c("No", "No", "Yes", "No", "Yes"),
         "District FE" = c("No", "Yes", "Yes", "No", "Yes")), 
       file = paste0(tab.out, "tableJ20.tex"))


# Women subset
controls_dn = c("age + sc + st ")
dn.lm1 = paste(base_dn) %>% as.formula() %>%
  lm_robust(data = dn.dta %>% filter(fem==1) %>%
              rename(vwreserv = reserved_fem, vstreserv = reserved_st),
            clusters = uniq_panch_code, se_type = "stata")

dn.lm2 = paste(base_dn) %>% as.formula() %>%
  lm_robust(data = dn.dta %>% filter(fem==1) %>%
              rename(vwreserv = reserved_fem, vstreserv = reserved_st),
            clusters = uniq_panch_code, se_type = "stata",
            fixed_effects = ~uniq_dist_code)

dn.lm3 = paste(base_dn, "+", controls_dn, "+", "prop_st") %>% as.formula() %>%
  lm_robust(data = dn.dta %>% filter(fem==1) %>%
              rename(vwreserv = reserved_fem, vstreserv = reserved_st),
            clusters = uniq_panch_code, se_type = "stata",
            fixed_effects = ~uniq_dist_code)

dn.dta = dn.dta %>%
  mutate(reserved_scst = ifelse(reserved_sc==1 | reserved_st==1, 1, 0))

dn.lm5 = paste(base_dn_scst) %>% as.formula() %>% lm_robust(data = dn.dta %>% filter(fem==1) %>%
                                                              rename(vwreserv = reserved_fem, vscstreserv = reserved_scst),
                                                            clusters = uniq_panch_code, se_type = "stata")

dn.lm6 = paste(base_dn_scst, "+", controls_dn) %>% as.formula() %>% lm_robust(data = dn.dta %>% filter(fem==1) %>%
                                                                                rename(vwreserv = reserved_fem, vscstreserv = reserved_scst),
                                                                              clusters = uniq_panch_code, se_type = "stata",
                                                                              fixed_effects = ~uniq_dist_code)




dv.means = rep(round(mean(dn.dta$marry_other_caste_bin[dn.dta$fem==1], na.rm = T), 3), 5)

screenreg(l = list(dn.lm1, dn.lm2, dn.lm3, dn.lm5, dn.lm6),
          include.ci = F, include.rmse = F,
          include.nobs = T,include.adjrs = FALSE, include.rsquared = F,
          digits = 3, include.nclust = F,
          omit.coef = c("age|^sc$|^st$|fem|Intercept|prop\\_st"),
          custom.coef.names = c("Women's Quota", "Caste Quota (ST)", "W X C (ST) Quota",
                                "Caste Quota (ST/SC)", "W X Caste (ST/SC) Quota"),
          stars = c(0.01, 0.05, 0.1),
          custom.model.names = c("Women", "Women", "Women", "Women", "Women"),
          custom.gof.rows = list(
            "Caste Quota Type" = c("ST", "ST", "ST", "ST/SC", "ST/SC"),
            "DV Mean" = dv.means, 
            "Controls" = c("No", "No", "Yes", "No", "Yes"),
            "District FE" = c("No", "Yes", "Yes", "No", "Yes")), 
          file = paste0(tab.out, "tableJ21.tex"))








# ST women subset
controls_dn = c("age ")


dn.lm1 = paste(base_dn) %>% as.formula() %>%
  lm_robust(data = dn.dta %>% filter(fem==1 & st==1) %>%
              rename(vwreserv = reserved_fem, vstreserv = reserved_st),
            clusters = uniq_panch_code, se_type = "stata")

dn.lm2 = paste(base_dn) %>% as.formula() %>%
  lm_robust(data = dn.dta %>% filter(fem==1 & st==1) %>%
              rename(vwreserv = reserved_fem, vstreserv = reserved_st),
            clusters = uniq_panch_code, se_type = "stata",
            fixed_effects = ~uniq_dist_code)

dn.lm3 = paste(base_dn, "+", controls_dn, "+", "prop_st") %>% as.formula() %>%
  lm_robust(data = dn.dta %>% filter(fem==1 & st==1) %>%
              rename(vwreserv = reserved_fem, vstreserv = reserved_st),
            clusters = uniq_panch_code, se_type = "stata",
            fixed_effects = ~uniq_dist_code)

dn.dta = dn.dta %>%
  mutate(reserved_scst = ifelse(reserved_sc==1 | reserved_st==1, 1, 0))

dn.lm5 = paste(base_dn_scst) %>% as.formula() %>% lm_robust(data = dn.dta %>% filter(fem==1 & (st==1 | sc==1)) %>%
                                                              rename(vwreserv = reserved_fem, vscstreserv = reserved_scst),
                                                            clusters = uniq_panch_code, se_type = "stata")

dn.lm6 = paste(base_dn_scst, "+", controls_dn) %>% as.formula() %>% lm_robust(data = dn.dta %>% filter(fem==1 & (st==1|sc==1)) %>%
                                                                                rename(vwreserv = reserved_fem, vscstreserv = reserved_scst),
                                                                              clusters = uniq_panch_code, se_type = "stata",
                                                                              fixed_effects = ~uniq_dist_code)

dv.means1 = rep(round(mean(dn.dta$marry_other_caste_bin[dn.dta$fem==1&dn.dta$st==1], na.rm = T), 3), 3)
dv.means2 = rep(round(mean(dn.dta$marry_other_caste_bin[dn.dta$fem==1& (dn.dta$st==1|dn.dta$sc==1)], na.rm = T), 3), 2)


screenreg(l = list(dn.lm1, dn.lm2, dn.lm3, dn.lm5, dn.lm6),
          include.ci = F, include.rmse = F,
          include.nobs = T,include.adjrs = FALSE, include.rsquared = F,
          digits = 3, include.nclust = F,
          omit.coef = c("age|^sc$|^st$|fem|Intercept|prop\\_st"),
          custom.coef.names = c("Women's Quota", "Caste Quota (ST)", "W X C (ST) Quota",
                                "Caste Quota (ST/SC)", "W X Caste (ST/SC) Quota"),
          stars = c(0.01, 0.05, 0.1),
          custom.model.names = c("ST Women", "ST Women", "ST Women", "ST/SC Women", "ST/SC Women"),
          custom.gof.rows = list(
            "Caste Quota Type" = c("ST", "ST", "ST", "ST/SC", "ST/SC"),
            "DV Mean" = c(dv.means1, dv.means2), 
            "Controls" = c("No", "No", "Yes", "No", "Yes"),
            "District FE" = c("No", "Yes", "Yes", "No", "Yes")), 
          file = paste0(tab.out, "tableJ22.tex"))


